Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.
library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed
library(bayesplot) # needed for MCMC diagnostics
library(DHARMa) # model validation
library(ggdist) # partial plots
library(tidybayes) # partial plots
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans
library(ggpubr) # arranging figuresThis data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty
rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |>
mutate(dev.temp = as.factor(dev.temp),
replicate = str_sub(sampleID, -1,-1),
population = factor(population)) |>
# number of observations = 5758
filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
sampleID != "60.LCKM.152.30.1" # 5637 - 76 = 5542
) |>
filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes)
group_by(sampleID) |>
mutate(max_value_index = which.max(rate.output),
row_number = row_number(),
max.rate = max(rate.output),
low.rate = mean(rate.output[order(rate.output)[1:10]]),
net.rate = max.rate - low.rate) |>
filter(row_number <= max_value_index) |>
ungroup() |>
mutate(cMASS = scale(mass, center=TRUE, scale=FALSE)[,1],
cAge = scale(age, center=TRUE, scale=FALSE)[,1],
region = factor(region),
dev.temp =factor(dev.temp),
LEVEL = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
tank >= 200 & tank <= 299 ~ 2,
tank >= 300 & tank <= 399 ~ 3,
TRUE ~ NA_real_)),
female = factor(female))eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figGreat lets start to fit our models. There are number of variables within our data frame some of which may or may not be important. To investigate which variables are important and which are not will test different hypotheses that use different combinations of variables that are suited to answer reasonable hypotheses.
From out exploratory data analysis we can see that due to the upwards infliction within the data as fish approach CTmax we will likely need to apply a second order polynomial function to out temperature covariate. Due to the large sample size we have we will also explore how will General Additive Models can be used to model our data.
We will start out by testing two different hypotheses.
model1 <- brm(rate.output2 ~ region*dev.temp*poly(Temperature, 2) + cMASS,
family=gaussian(),
data=lat_resp_dat2,
warmup = 4000,
iter=10000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
chains=4,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## No divergences to plot.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model.equip <- brm(rate.output2 ~ region*dev.temp*poly(Temperature, 2) + cMASS + system + chamber,
family=gaussian(),
data=lat_resp_dat2,
warmup = 4000,
iter=10000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
chains=4,
thin = 5,
control = list(adapt_delta=0.95))## Compiling Stan program...
## Start sampling
## No divergences to plot.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## Warning:
## 2 (0.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
## Warning:
## 2 (0.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
## elpd_diff se_diff
## model.equip 0.0 0.0
## model1 -54.1 10.6
## Output of model 'model1':
##
## Computed from 4800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 623.0 63.1
## p_loo 25.6 1.4
## looic -1246.0 126.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model.equip':
##
## Computed from 4800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo 677.1 64.4
## p_loo 30.0 1.6
## looic -1354.2 128.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.1]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model.equip 0.0 0.0
## model1 -54.1 10.6
It seems like additional information of chamber and system do not produce a more informative model
Before we validate the model we will include out random effects:
findpriors = get_prior(rate.output2 ~ region*dev.temp*poly(Temperature, 2) + cMASS + (1|LEVEL) + (1|female),
family=gaussian(),
data=lat_resp_dat2); findpriorspriors = c(set_prior("normal(0,100)", class="b", coef="cMASS"),
set_prior("normal(0,100)", class="b", coef="dev.temp30"),
set_prior("normal(0,100)", class="b", coef="dev.temp30:polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="dev.temp30:polyTemperature22"),
set_prior("normal(0,100)", class="b", coef="dev.temp31"),
set_prior("normal(0,100)", class="b", coef="dev.temp31:polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="dev.temp31:polyTemperature22"),
set_prior("normal(0,100)", class="b", coef="polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="polyTemperature22"),
set_prior("normal(0,100)", class="b", coef="regionleading"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp30"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp30:polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp30:polyTemperature22"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp31"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp31:polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="regionleading:dev.temp31:polyTemperature22"),
set_prior("normal(0,100)", class="b", coef="regionleading:polyTemperature21"),
set_prior("normal(0,100)", class="b", coef="regionleading:polyTemperature22")
)
model1.re <- brm(rate.output2 ~ region*dev.temp*poly(Temperature, 2) + cMASS + (1|LEVEL) + (1|female),
family=gaussian(),
prior = priors,
data=lat_resp_dat2,
warmup = 1500,
iter=4000,
seed=123,
cores=4,
save_pars = save_pars(all=TRUE),
chains=3,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 31 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
We will continue with the first model. lets start validating the model to see how it performed. The initial diagnostics seemed to be pretty good. Now lets check to see if it is violating any assumptions.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
preds <- posterior_predict(model1.re, ndraws=250, summary=FALSE)
model1.re.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = na.omit(lat_resp_dat2$rate.output2),
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = 'student')
plot(model1.re.resids) ## Warning: There were 31 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output2 ~ region * dev.temp * poly(Temperature, 2) + cMASS + (1 | LEVEL) + (1 | female)
## Data: lat_resp_dat2 (Number of observations: 4776)
## Draws: 3 chains, each with iter = 4000; warmup = 1500; thin = 5;
## total post-warmup draws = 1500
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.13 0.03 0.09 0.20 1.00 1261 1316
##
## ~LEVEL (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.07 0.01 0.26 1.00 1194 906
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 0.62 0.05 0.51 0.73
## regionleading 0.03 0.07 -0.12 0.17
## dev.temp30 0.03 0.01 0.01 0.05
## dev.temp31 0.06 0.01 0.04 0.08
## polyTemperature21 8.04 0.47 7.12 8.96
## polyTemperature22 5.03 0.49 4.08 6.00
## cMASS 416.85 6.01 405.31 429.39
## regionleading:dev.temp30 -0.06 0.01 -0.08 -0.03
## regionleading:dev.temp31 0.00 0.01 -0.02 0.03
## regionleading:polyTemperature21 -0.57 0.69 -1.96 0.79
## regionleading:polyTemperature22 0.77 0.70 -0.61 2.13
## dev.temp30:polyTemperature21 -0.18 0.65 -1.49 1.08
## dev.temp31:polyTemperature21 -1.07 0.63 -2.26 0.18
## dev.temp30:polyTemperature22 -0.65 0.64 -1.90 0.57
## dev.temp31:polyTemperature22 -1.10 0.65 -2.34 0.17
## regionleading:dev.temp30:polyTemperature21 0.45 1.00 -1.43 2.41
## regionleading:dev.temp31:polyTemperature21 0.15 0.93 -1.70 1.96
## regionleading:dev.temp30:polyTemperature22 1.19 0.99 -0.68 3.20
## regionleading:dev.temp31:polyTemperature22 0.50 0.93 -1.30 2.32
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1254 1076
## regionleading 1.00 1376 1417
## dev.temp30 1.00 1359 1502
## dev.temp31 1.00 1453 1420
## polyTemperature21 1.00 1556 1426
## polyTemperature22 1.00 1341 1383
## cMASS 1.00 1418 1240
## regionleading:dev.temp30 1.00 1377 1520
## regionleading:dev.temp31 1.00 1257 1422
## regionleading:polyTemperature21 1.00 1456 1457
## regionleading:polyTemperature22 1.00 1280 1419
## dev.temp30:polyTemperature21 1.00 1483 1268
## dev.temp31:polyTemperature21 1.00 1402 1407
## dev.temp30:polyTemperature22 1.00 1385 1224
## dev.temp31:polyTemperature22 1.00 1470 1499
## regionleading:dev.temp30:polyTemperature21 1.00 1427 1482
## regionleading:dev.temp31:polyTemperature21 1.00 1523 1478
## regionleading:dev.temp30:polyTemperature22 1.00 1468 1367
## regionleading:dev.temp31:polyTemperature22 1.00 1457 1342
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.19 0.00 0.19 0.19 1.00 1479 1443
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#model1.re.wo |> get_variables()
model1.re |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()